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INTRODUCTION 


Application of three-dimensional design and analysis codes to aircraft design 
problems has become increasingly important. Reduced development time and costs are 
the most obvious benefits. Also, these new computational methods may make it possi- 
ble to design configurations with improved performance, particularly in the transonic 
speed range. Since three-dimensional transonic inverse codes have not been fully 
developed, the trial-and-error design process can still be quite lengthy even with 
good analysis codes. Optimization can be applied to shorten the design cycle. 

The purpose of this study is to utilize a three-dimensional transonic analysis 
method coupled with a numerical optimization procedure to reduce the wave drag of a 
fighter wing at transonic maneuver conditions. The code used in this study is called 
PANDORA ^Preliminary Automated *tonierical Design of Realistic JJarcraft). Developed by 
Aidala (ref. 1), the code combines a three-dimensional transonic analysis method and 
a numerical optimization procedure. The analysis method is a modified version of 
Boppe's small disturbance transonic wing-body code (ref. 2), and the optimization 
procedure is that of Vanderplaats (ref. 3). Wave drag is evaluated through the use 
of a formula based on the loss in momentum across an isentropic shock. Ihe optimiza- 
tion procedure minimizes wave drag by modifying the wing section contours from root 
to tip via a wing profile shape function. Wing angle of attack is also allowed to 
vary during the optimization procedure. 

An existing fighter wing was redesigned with the PANDORA code to provide a sig- 
nificant reduction in wave drag at transonic maneuver conditions. This paper dis- 
cusses this redesign effort and presents results from optimization. Comparisons 
between the trends in the airfoil modifications produced by the optimization proce- 
dure and those that resulted previously from applying various transonic analysis 
codes are also presented. 


SYMBOLS 

AR aspect ratio 

b wing span 

C n • wing induced drag coefficient 

C D p wing pressure drag coefficient 

C L total configuration lift coefficient 

Cp pressure coefficient 

c local chord 

c av wing average chord 

c d p wing section pressure drag coefficient 

E Oswald efficiency factor 


p local pressure 

t/2c half-thickness ratio 

u local streamwise velocity component 

x streamwise direction 

x/c chord fraction 

y spanwise direction 

z vertical direction 

z/c camber line 

Az/c change in camber line 

ri wing span station, 2y/b 

p density 

<|> perturbation velocity potential 


DISCUSSION 

Transonic Analysis Method 

The PANDORA code used in this study was developed by Aidala (ref, 1) and con- 
sists of a three-dimensional transonic analysis method coupled with a numerical opti- 
mization procedure. The analysis method is essentially Boppe's wing-body code 
(ref. 2) modified to include the capability of modeling a canard. The governing 
potential flow equation used in the analysis code is in an extended small-disturbance 
form. Extra terms have been added to improve the resolution of swept shock waves and 
to obtain a better approximation of the critical velocity. The flow solution is 
obtained in a multiple embedded grid system using successive line overrelaxation. 

The grid scheme and solution process are described briefly. References 1 and 2 pro- 
vide more detailed information. 

The computational space used in the embedded grid scheme consists of a crude 
Cartesian mesh and embedded fine grids (fig. 1). The physical domain whose bound- 
aries correspond to infinity is transformed to the computational domain with finite 
boundaries. At the bounding planes, the potential is set to zero except at the down- 
stream and symmetry planes. The flow-field potential is obtained at the downstream 
plane from the equation 


<b +6 = 0 
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At the symmetry plane, the conditions 


4> = o 
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<b =0 
xy 

are applied* 

A secondary mesh system consisting of individual fine grid arrays for the wing 
and canard is embedded in the crude grid. The fine grid system is sheared and 
tapered to fit the wing and canard planforms. Wing and canard fine grid arrays are 
located where crude streamwise mesh planes intersect the planforms. These fine grids 
are evenly spaced in the vertical and spanwise directions. 

The crude and the fine grid solutions interact by alternately updating the solu- 
tion on the crude and fine grids. Initially, calculations are made only on the crude 
grid in order to develop the global characteristics of the flow field. The interac- 
tion scheme then proceeds by updating the potentials on the perimeter of the fine 
grids with potentials determined by the previous crude grid solution. Similarly, the 
crude grid potentials on the lifting surface are determined by the previous fine grid 
solutions. This interaction process continues until the solution has converged 
sufficiently. 

Output from the analysis code consists of pressure coefficients, convergence 
history, spanwise load, moment and drag distributions, and total configuration force 
and moment coefficients. 


Optimization Procedure 

Through a series of subroutines, the analysis method is coupled with the optimi- 
zation procedure CONMIN developed by Vanderplaats (ref. 3). This procedure is struc- 
tured to solve linear and nonlinear constrained minimization problems by a modified 
method of feasible directions. The design task presented in this paper serves as an 
example demonstrating the optimization process. Reference 3 provides more detailed 
information • 

The current design task attempts to minimize wave drag at a specific lift coef- 
ficient. In optimization terminology, wave drag is the objective function, and the 
lift coefficient is the lower bound constraint. A set of independent design vari- 
ables are defined such that they control the shape of the wing surface contour and 
the wing angle of attack. The optimization process is provided with reference values 
for the objective function and constraint by the analysis method. These reference 
values are used throughout the optimization process as discussed below. 

The optimization process begins with the sequential perturbation of each design 
variable by 10 percent to determine the effects on wave drag (objective function) and 
lift (constrained variable). These effects are based on the changes in the objective 
function and constrained variable relative to their reference values. From this 
information, the gradients of both the objective function and constrained variable 
are calculated with respect to the design variables. With these gradients, a search 
direction is determined, based on superposition, which will reduce the objective 
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function if there is no adverse interaction between the variables. After a search 
direction is determined, the design variables are incremented (stepped) accordingly. 
For clarity, this portion of the optimization procedure will be referred to as 
step 1 . 

In step 2, the process returns to the analysis code with the incremented design 
variables to calculate the resulting values of the objective function and constrained 
variable. A comparison is then made between these new values and the reference val- 
ues (step 3). If the objective function increased or the constraint has been vio- 
lated, the step size in the search direction is decreased and the process returns to 
step 2. If the objective function decreased without violating the constraint, the 
design variables are incremented, and steps 2 and 3 are repeated. When the objective 
function can no longer be reduced without violating the constraint, optimization is 
terminated, and a local optimum is assumed to have been found. This constitutes one 
optimization iteration. The user may specify more optimization iterations, in which 
case the gradients would be recalculated to determine a new search direction. 

For the wave-drag-minimization problem presented in this paper, only one optimi- 
zation iteration was desired; that is, only one search direction was investigated. 

It is fundamental to this optimization procedure that a converged analysis solu- 
tion be obtained prior to initiating optimization. Since the solution for the ini- 
tial geometry and freestream conditions is the starting solution for all analysis 
computations performed during optimization, the comparisons described in steps 1 
and 3 above are valid if the starting solution is converged. Without this require- 
ment, the values of the objective function and constrained variable produced with 
incremented design variables could be from a more converged solution than their ref- 
erence values • 

During the study, a problem of slow convergence of the analysis solution arose 
as a result of high loading at the wing tip. To overcome this problem, a separate 
analysis run was used to determine the reference values of the objective function and 
constraint. The number of iterations for the reference analysis run was determined 
by the total number of iterations used to develop the solution for the perturbed 
geometry. As a result, valid comparisons were made between the initial geometry and 
the perturbed geometry without satisfying the convergence requirement. It should be 
noted that initially it was believed that only the reference value for the objective 
function needed to be determined in this manner. This procedure also permitted com- 
parisons between the initial and optimized wing contours at the same time step in the 
convergence history. 


Wave Drag Computation 

The wave drag coefficient is calculated in the original version of PANDORA by 
subtraction of the induced drag coefficient from the wing drag coefficient. The wing 
drag coefficient is calculated by integration of the section drag coefficient along 
the span 
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and the induced drag coefficient is obtained from the equation 


C 


D,i 



TT (AR)E 


In place of this determination of wave drag, a more direct computational method 
was sought that could be used as the objective function for optimization. Steger and 
Baldwin (ref. 4) have shown that a realistic wave drag quantity can be computed by 
relating wave drag to the loss in momentum across an isentropic shock. Prom the 
relation 

P 2 + P 2 U 2 2 < P 1 + P 1 U 1 2 


where subscript 2 denotes conditions behind a normal shock 
conditions in front of a normal shock, an equation for the 
area at a given point on the shock can be derived in terms 


and subscript 1 denotes 
local wave drag per unit 
of the loss in momentum as 


Ad = (p 
wave 1 


D u 
M 1 1 


•> - <p 2 + 


P 2 U 2 


Note that no account is taken of the shock sweep. Integration of this equation gives 
local wave drag at each wing span station. The total wave drag coefficient is then 
obtained through integration and nondimensionalization of the local values across the 
span. 


Implementation of this wave drag computation requires a reliable method to 
determine the shock location at each wing span station. Two methods were 
investigated. 

The first method determined shock location based on the local flow velocities. 
The location of the first subsonic point aft of the supersonic zone was used as the 
shock location. This method has been shown to work for small disturbance treatment 
of two-dimensional transonic flow in reference 4. However, from figure 2 (method 1), 
it can be seen that for the present three-dimensional problem there is a significant 
difference between the onset of the shock and the shock location as determined by 
this method. These results indicate that the first method is not a reliable indi- 
cator of shock location for the present problem. 

In the second method, pressure coefficients were calculated at each span station 
starting from the leading edge and proceeding in the streamwise direction to the 
first subsonic point aft of the supersonic zone. The gradients of the local pressure 
coefficients with respect to chord location were then computed. The chord location 
of the largest pressure gradient was deemed the shock location. Figure 2 (method 2) 
illustrates the computed shock location obtained with the second method for several 
span stations. To initiate the wave drag computation at the onset of the shock, it 
was necessary to begin calculations four mesh spaces (0.04c) upstream of this point. 
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This was a sufficient amount , since the method consistently selected the second or 
third mesh point on the shock regardless of shock strength (fig. 2). This method of 
determining shock location was used in the wave drag procedure. 


Wing Geometry /Shape Function 

The initial geometry is an early version of the Pathfinder II wing-body fighter 
configuration (fig. 3) developed at the Langley Research Center (ref. 5). The wing 
planform has a leading-edge sweep of 45° and a trailing-edge sweep of 11.9°. Aspect 
ratio is 3.28, and taper ratio is 0.214. An axi symmetric model of the fuselage is 
used. A second wing design developed from the initial geometry through the use of 
various transonic analysis codes was used to evaluate the trends of the optimized 
wing profile shapes. Figure 4 gives the thickness distribution for both wing geome- 
tries. The camber lines for the initial wing geometry are given in figure 5, and the 
differences in camber lines between these two configurations are given in figure 6. 

It was anticipated that through optimization the wing profile shapes of the initial 
geometry would be modified such that the trends observed in figure 6 would be 
produced. 

A shape function was constructed such that the entire wing could be optimized in 
one optimization run. To make the resulting wing more feasible, six variables which 
allow emphasis on specific areas of the wing where modifications are more important 
were included in this shape function. The shape function is 

2 

■ 4 [ C 2 V» a , + C 3 V" 1 a 2 + C 4 V"’ a 3 ](t - *o) 
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where a value of 0.4 was taken for x Q , and Az represents the change in the camber 
line. No changes were made to the wing profile shapes forward of 0.4c. The thick- 
ness distribution was held fixed so it was not necessary to optimize the lower sur- 
face separately. The variables C 2 to are design variables whose values are 

determined by the optimization procedure. The coefficients A^ and represent 

wing changes which are largest at the root and decrease to zero at the tip; the coef- 
ficients A^ and B 2 represent wing changes which are uniform from root to tip; and 
the coefficients A^ and B^ represent wing changes which are largest at the tip 
and decrease to zero at the root. (See table I.) Coefficients A^ to A^ and 
B^ to B 3 may be "turned on or off" by specifying each coefficient a^ to 
and b 1 to b^ to be either one or zero, to control the type of wing change 
allowable. 

Two sets of values for the a and b coefficients were constructed for optimi- 
zation. (See table II.) For variable set 1, only a 2 and b ^ were assigned values 
of one so that a uniform change would be made across the span with slightly more 
emphasis on the root sections. Because the remaining a and b coefficients were 
set to zero, only two design variables were used in the shape function. The effect 
of variable set 2 was to gradually increase the wing modifications from the root to 
the tip. Four design variables were used in this shape function. 
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RESULTS 


The starting conditions for the optimization procedure were an angle of attack 
of 12.6° and a freestream Mach number of 0.85 at a lift coefficient of 0.562. The 
lower bound constraint on lift was set to this value. As discussed previously, a 
reference analysis run was used to determine the reference value for wave drag to 
ensure that comparisons were made at the same time step in the convergence history. 
The reference value for wave drag coefficient was determined to be 0.0077. One opti- 
mization iteration was attempted with 60 crude/fine iterations allowed for each flow 
solution requested by the optimization procedure. Variable set 1, given in table II, 
was used for the wing shape function, so that two design variables for wing shape 
modifications and one design variable for changes in wing angle of attack were 
employed. This three-design-variable optimization produced a lift coefficient of 
0.561 and a wave drag coefficient of 0.0070. The wing angle of attack increased to 
12.7°. The optimized wing did lower the wave drag coefficient but not without reduc- 
ing the lift coefficient. Virtually identical modifications were made across the 
span because of design variable C^, which was more influential than in the 

shape function. The resulting modifications made to the initial geometry, shown in 
figure 7, did not reflect the trend that was anticipated. This indicated that more 
emphasis should be placed on modifying the outboard stations. It also appeared that 
more than two design variables should be used for the wing shape modifications. 
Evaluation of the gradients of the objective function and constraint made it apparent 
that a reference value for the lift coefficient should be determined and used in the 
same manner as the reference value for the wave drag. Because the lift coefficient 
had not sufficiently converged prior to this optimization attempt, the efficiency of 
the search direction may have been reduced. 

A second optimization was performed starting from a different set of initial 
conditions. These conditions were changed to reduce the loading at the wing tip 
while maintaining a higher lift coefficient. These conditions were an angle of 
attack of 16° at a freestream Mach number of 0.80. A reference analysis run gave the 
reference values of 0.0077 for wave drag coefficient and 0.716 for lift coefficient. 
Variable set 2, given in table II, was constructed such that emphasis would be placed 
on modifying outboard wing stations and increasing the number of design variables 
influencing the shape function from two to four. One optimization iteration of the 
configuration yielded a wave drag coefficient of 0.0081 and a lift coefficient of 
0.755. The wing angle of attack was reduced to 15.6°. Figure 8 illustrates the 
modifications to the wing profiles which reflect the desired trend. Although wave 
drag increased, there was also an increase in the lift coefficient. Without further 
modifications to the wing profiles, wave drag could be reduced by lowering the wing 
angle of attack until the lift coefficient reaches its constraint value. To accom- 
plish wing profile modifications and additional reduction in wing angle of attack 
through optimization would have required more than one optimization iteration; that 
is, more than one search direction would have to be interrogated. 

Through analysis runs, the lift coefficients of the optimized and initial geom- 
etry were matched to determine the optimized angle of attack. To obtain more con- 
verged solutions, these analysis runs were allowed to run longer than those pre- 
viously used. After a total of 560 iterations at the same conditions as the second 
optimization case, the lift coefficient of the initial geometry was 0.735, and the 
wave drag coefficient was 0.0100. A lift coefficient of 0.738 for the optimized 
geometry was obtained by reducing the wing angle of attack to 13°. This resulted in 
a wave drag coefficient of 0.0040. By matching the lift coefficients, a significant 
decrease in wave drag was achieved. A comparison of the pressure distributions for 
the baseline and optimized geometries at the above conditions is shown in figure 9. 
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The shock strength and leading-edge pressure peaks were significantly reduced across 
the span. Except in the tip region, the onset of the shock moved aft approximately 
0 • 08c . 


CONCLUDING REMARKS 

The original computational method for determining wave drag in a three- 
dimensional transonic analysis method was replaced by a wave drag formula based on 
the loss in momentum across an isentropic shock. This formula was used as the objec- 
tive function in an optimization procedure coupled with the analysis method to reduce 
the wave drag of a fighter wing at transonic maneuver conditions. The optimization 
procedure minimized wave drag at a specific lift coefficient by modifying the wing 
section contours defined by a wing profile shape function. Through the use of the 
three-dimensional transonic analysis code and optimization procedure, a fighter wing 
was redesigned with a significant reduction in wave drag at transonic maneuver condi- 
tions. Leading-edge pressure peaks and shock strength were significantly reduced 
while maintaining a high lift coefficient. 

Although a successful optimization was performed with wave drag as the objective 
function, the high sensitivity to solution convergence demonstrated by this computa- 
tion may have reduced the efficiency of the optimization process. Further study is 
needed to reduce the effects of convergence on the wave drag computation. In addi- 
tion, shock sweep effects should be included in the computation. 

The shape function used by the optimization procedure to modify the wing section 
contours enabled the optimization procedure to modify all wing section contours in 
one optimization iteration. However, having this capability does put severe limita- 
tions on the optimization process. For the particular case presented here, the types 
of wing contour modifications necessary to reduce wave drag were produced through the 
use of this shape function. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
December 21, 1983 
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